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Abstract: This paper presents a new event generator, ALPGEN, dedicated to the study 
of multiparton hard processes in hadronic collisions. The code performs, at the leading 
order in QCD and EW interactions, the calculation of the exact matrix elements for a large 
set of parton-level processes of interest in the study of the Tevatron and LHC data. The 
current version of the code describes the following final states: (W — ► ff')QQ + N jets (Q 
being a heavy quark, and / = £,q), with N < 4; (Z/j* — >■ ff)QQ + N jets (/ = £,v), 
with N < 4; (W -> //') + charm + N jets (/ = £,q, N < 5); (W -» //') + N jets 
(/ = £, q) and (Z/j* -> ff) + N jets (/ = £, u), with N < 6; nW + mZ + IH + N jets, 
with n + m + l + N<8, N < 3, including all 2-fermion decay modes of W and Z bosons, 
with spin correlations; QQ + N jets, with t — > bff decays and relative spin correlations 
included if Q = t, and iV < 6; QQQ'Q + N jets, with Q and Q' heavy quarks (possibly 
equal) and N < 4; HQQ + N jets, with t — ► bff decays and relative spin correlations 
included if Q = t and N < 4; N jets, with N < 6. Parton-level events are generated, 
providing full information on their colour and flavour structure, enabling the evolution of 
the partons into fully hadronised final states. 
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1. Introduction 



The large energies available in current and forthcoming hadronic colliders make final states 
with several hard and well separated jets a rather common phenomenon. These multijet 
final states can originate directly from hard QCD radiative processes [1], or from the decay 
of massive particles, such as for example W and Z gauge bosons, top quarks, Higgs bosons, 
supersymmetric particles, etc. Whether our interest is in accurate measurements of top 
quarks or in the search for more exotic states [2], multijet final states always provide an 
important observable, and the study of the backgrounds due to QCD is an essential part 
of any experimental analysis. 

Several examples of calculation of multijet cross-sections in hadronic collisions exist 
in the literature. Some of them are included in parton-level Monte Carlo (MC) event 
generators, where final states consisting of hard and well isolated partons are generated. 
Among the most used and best documented examples are PAPAGENO [3] (a compilation of 
several partonic processes), VECBOS [4] (for production of W/Z bosons in association with 
up to 4 jets), NJETS [5] (for production of up to 6 jets). The range of jet multiplicities 
calculable in practice for purely QCD processes in hadronic collisions has recently been 
extended in [6], where new techniques [7] to deal with the complexity of the multijet flavour 
configurations have been implemented in the calculation of up to 7 jets. Finally, programs 
for the automatic generation of user-specified parton-level processes exist and have been 
used in the literature for the calculation of many important reactions in hadronic collisions: 
MADGRAPH [8], CompHEP [9], GRACE [10] and AMEGIC++ [11]. 

In order to use these results in practical analyses of the experimental data, the cal- 
culations need to be completed with the treatment of the higher-order corrections leading 
to the development of partonic cascades, and with the subsequent transformation of the 
partons into observable hadrons. MC programs such as HERWIG [12], PYTHIA [13] or 
ISAJET [14] are available to carry out these last two steps. The consistent combination of 
the parton-level calculations with the partonic evolution given by the shower MC programs 
is the subject of extensive work. Several approaches to this problem have been proposed in 
the case of low-order processes, where one is interested in final states with at most one ex- 
tra jet relative to a given Born-level configuration (for example W, Z plus jet [15] and ti + 
jet [16]). In these cases, the proposed algorithms correct the probability for hard-emission 
estimated by the approximation of the shower-evolution programs, using the value of the 
exact real-emission higher-order matrix element. In other studies [17, 18, 19, 20], cover- 
ing at this time jet emission in association with DY, vector boson pairs and heavy quark 
pairs, the full set of virtual and real NLO corrections to the partonic matrix elements has 
been merged with the HERWIG MC. See in particular [18] for a complete discussion of the 
problem of matching NLO parton level and leading-logarithmic (LL) shower generators. 
Some of the problems raised by a consistent matching of NLO parton level calculations 
and next-to-leading logarithm (NLL) shower evolution are discussed in [21]. 

In the case of large jet multiplicities, the complexity of the matrix element evaluation 
and of its singularity structure prevents so far the application of these approaches. Recently, 
a new strategy has been introduced [22, 23] for the merging of multijet matrix elements 
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with the shower development. This involves a reweighting of the matrix element weights 
with Sudakov form factors, and the veto of shower emissions in regions of phase-space 
already covered by the parton-level configurations. After the shower evolution, samples 
of different parton-level multiplicity are combined together to obtain inclusive samples of 
arbitrary jet multiplicity, double counting being limited to subleading effects. A complete 
application of these ideas has been achieved for inclusive hadronic final states in e + e~ 
collisions [24, 23], and a proposal has been formulated for the extension to the hadron 
collider case [25]. Explicit implementations are being developed [26]. 

We discussed in [27] a theoretical framework for the evaluation and generation of 
events in such a way as to enable the subsequent perturbative evolution using a shower 
MC program. In [28] we presented a complete application to the case where a W boson 
is produced in association with a heavy quark pair, plus up to 4 additional light partons. 
The relative MC code allows a complete description of these final states, from the leading 
order (LO) matrix element computation to the perturbative evolution and hadronization 
carried out using HERWIG. 

We recently extended the application of the ideas contained in [27], completing a library 
of MC codes for hadronic collisions including the following new processes: WQQ + N jets 
and Z/'y* QQ + N jets (Q being a heavy quark), with N < 4; W + charm + N jets; 
W + N jets and Z + N jets (N < 6); nW + mZ + IH + N jets, with n + m + / + A<8and 
N < 3; QQ + N jets (N < 6); QQQ'Q' + N jets, with Q and Q' heavy quarks (possibly 
equal) and N < 4; HQQ + N jets (N < 4); N jets, with N < 6. Details on the decay mode 
options for the various unstable particles will be given below. For all of these processes, an 
interface to both HERWIG and PYTHIA is provided. In the present work, we document the 
contents and use of this library. The emphasis is on the code itself and not on the physics 
approach, which is discussed in more detail in [28]. Numerical results for some benchmark 
processes are nevertheless presented. 

Independent work on the merging of parton-level calculations with shower MC's has 
been pursued by the CompHEP group [29], by the GRACE group [30] and, more recently, 
in [31, 32, 33]. To the best of our knowledge, a large fraction of the matrix element 
calculations documented in this paper have however never been performed before. 

Section 2 reviews the general structure of the codes, covering both the aspects of the 
parton-level calculations, and of the shower evolution. Section 3 discusses the features of 
the implementation of each individual hard process in our library. A technical Appendix 
will provide more explicit details on the programs and their use. 

The library containing all codes described in this work can be downloaded from the 
following URL: http : //mlm . home . cern . ch/mlm/ alpgen . 

2. The general structure of the program 

The program consists of several building blocks (see the Appendix for details). A section of 
the code library defines the overall infrastructure of the generator, implementing the logical 
sequence of operations in a standard set of subroutine calls; this part is independent of the 
hard process selected and, among other things, it includes the algorithms needed for the 
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evaluation of the matrix elements, for the evaluation of the parton densities and for the 
bookkeeping and histogramming of the results. Each hard process has viceversa a separate 
set of code elements, which are specific to it. These include the process initialization, the 
phase-space generation, the extraction of flavour and colour structure of the event, and the 
default analysis routines. Each hard process corresponds to a specific executable, obtained 
by linking the relative process-dependent code elements with the process-independent ones. 
In addition to the above, a section of the code library deals with the shower evolution. As 
explained below, the shower evolution is performed as an independent step, following the 
generation of a sample of unweighted parton-level events. The code elements relative to this 
phase of the computation include the HERWIG and PYTHIA codes and the algorithm needed 
to transform the partonic input into a format which can be interpreted and processed by 
HERWIG and PYTHIA. We adopt the formatting standard proposed in the so called Les 
Houches accord [34]. 

As alluded to above, the program has two main modes of operation. In the first mode 
the code performs the parton-level calculation of the matrix elements relative to the selected 
hard processes, generating weighted events. Each weighted event is analysed on-line in a 
routine where the kinematics of the event can be studied, and histograms filled. The user 
has direct access to this analysis routine, and can adapt it to his needs. At the end of the 
run, differential distributions are obtained. A histogramming package is included in the 
code library, which generates a graphic output in the form of topdrawer [35] plots. This 
mode of running can also be used to easily get total cross sections in presence of some 
overall generation cuts (e.g. rates for production of jets above a given threshold), without 
looking at any particular differential distribution. 

In the second mode of operation the code generates unweighted parton level events and 
stores them to a file, for subsequent evolution via the parton-shower part of the program. 
The generation of parton-level events, and their shower evolution, are performed in two 
different phases by different code elements. Since the generation of unweighted parton-level 
events is typically by far the most CPU intensive component of the calculation, the storage 
of unweighted events allows to build up event data sets which can then be used efficiently 
for studies of hadronization systematics or realistic detector simulations. In this mode of 
operation the matrix-element calculation generates all the flavour and colour information 
necessary for the complete shower evolution. The kinematical, flavour and colour data for 
a given event are stored in a file, and are read in by the shower MC, which will process the 
event. In the rest of this Section we present in more detail these two running modes of the 
code. 

2.1 Parton-level generation and cross section evaluation 

In a nutshell, the calculation of the cross section for a given hard process is performed in 
the following steps: 

• The parameters required to define the hard process are passed to the code. These 
include the selection of jet multiplicity, the mass of possible heavy quarks, rapidity 
and transverse momentum cuts, etc. 
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• A first set of phase-space integration cycles is performed, with the goal of exploring 
how the cross-section is distributed in phase-space and among the possible contribut- 
ing subprocesses. Event by event, the following steps take place: 

- one subprocess (see later) is randomly selected; 

- a point in phase-space is randomly selected, consistent with the required kine- 
matical acceptance cuts; 

- the initial-state parton luminosity is evaluated for the chosen subprocess, and 
one among the possible flavour configurations is selected (see later) ; 

- spin and colour for each parton are randomly assigned; 

- the matrix element is evaluated, and the weight of the event is obtained after 
inclusion of the phase-space and parton-luminosity factors. A bookkeeping of 
the weights is kept for each individual subprocess and phase-space subvolume. 

• At the end of the first integration iteration, a map of the cross-section distribution 
among the different subprocesses and in phase-space is available. It will be used for 
subsequent integration cycles, where the phase-space and subprocess random sam- 
pling will be weighted by the respective probability distributions. 

• After the completion of a series of warm-up integration cycles (whose number is 
specified at the beginning of the run by the user) , the optimised integration grids are 
stored in a file. A final large-statistics run is then performed. After the generation 
of each event, its kinematics is analysed and histograms are filled. 

We shall now discuss in more detail the individual steps outlined above and the ingredients 
of the calculation. 

2.1.1 Selection of the subprocess 

The calculation of the cross section for multiparton final states involves typically the sum 
over a large set of subprocesses and flavour configurations. For example, in the case of 
WQQ + 2 jets we have, among others, the following subprocesses: 

qq' - WQQgg, qg - WQQgq', gq -v WQQgq', gg - WQQqq', qq' - WQQq"f, etc. 

(2.1) 

Each of these subprocesses receives contributions from several possible flavour configura- 
tions (e.g. ud — ► WQQgg , us — ► WQQgg, etc.). Our subdivision in subprocesses is 
designed to allow to sum the contribution of different flavour configurations by simply 
adding trivial factors such as parton densities or CKM factors, which factorize out of a 
single, flavour independent, matrix element. For example the overall contribution from the 
first process in the above list is given by 

\u1d~2 cos 2 6 C + u±S2 sin 2 9 C + c\s~2 cos 2 C + c\d,2 sin 2 6 C ] x \M(qq' -► WQQgg)] 2 , (2.2) 
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where qi = f(xi), i = 1,2, are the parton densities for the quark flavour q. Contributions 
from charge-conjugate or isospin-rotated states can also be summed up, after trivial mo- 
mentum exchanges. For example, the same matrix element calculation is used to describe 
the four events: 

u{pi)d(p 2 ) -> b(p 3 )b(p 4 )g(p 5 )g(p 6 )e + (p 5 )i/(p 6 ) 
u(pi)d(p 2 ) -> b(pz)b(p A )g(ps)g(ps)e-(p b )v(ps) 
d(pi)u{p 2 ) -> b{p'i)b{pA)g{p^)g{p( i )u{p5)e + (p 6 ) 
d{p\)u(p 2 ) -> b(p 3 )b(p 4 )g(p 5 )g(p 6 )v(p 5 )e- (p 6 ) . 

Event by event, the flavour configuration for the assigned subprocess is then selected with a 
probability proportional to the relative size of the individual contributions to the luminosity, 
weighted by the Cabibbo angles. 

Typically, only few among all possible subprocesses give a substantial contribution to 
the cross section. For each event, instead of summing the weight of all subprocesses, we 
calculate only one. This is selected with uniform probability during the first integration cy- 
cle, and a record is kept of each individual contribution to the cross-section. In subsequent 
integration iterations, the accumulated rates of the single channels are used to weight their 
selection probabilities. This will significantly improve the CPU performance of the code, 
and the unweighting efficiency. 

Due to the rapid growth in the number of subprocesses when quarks are added [7], 
we limited ourselves to processes with at most 2 pairs of light quarks (plus pairs of heavy 
quarks, when required). In all of the cases considered this is not, however, a significant 
limitation to the accuracy of the results. The full list of subprocesses available for each 
hard process is given in Section 3. 

2.1.2 Phase-space sampling 

The phase-space generation is optimised for each individual hard process, using generation 
variables which are most suitable to the application of typical hadron collider selection 
cuts. The phase-space is mapped with a multidimensional grid, and during the integration 
a record is kept of the total weight accumulated within each bin of the grid. To further 
contribute to the efficiency of the phase-space sampling, independent grids are employed 
to sample different subprocesses. In particular, we shall associate one phase-space grid to 
each of the following initial states: 

1. qq, qq' and charge conjugates 

2. qg and qg 

3. gq and gq 

4. gg 

5. qq, qq' and charge conjugates. 
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For example, the processes qq' 
grid. 



WQQgg and qq' — ► WQQqq share the same phase-space 



2.1.3 Matrix element calculation 

The calculation of the LO matrix elements for the selected hard process is performed 
using the ALPHA[36] algorithm, extended to QCD interactions as described in [27, 28]. As 
explained in detail in [28], use of the ALPHA algorithm is in our view essential in order 
to cope with the complexity of the problem. All mass effects are included in the case of 
massive quarks. The calculations are done after having selected, on an event-by-event basis, 
polarization, flavour and colour configurations, in order to be able to provide the shower 
MC's with the information necessary for the shower evolution. The sum over polarizations 
and colours is performed by summing over multiple events, in a MC fashion. The choice 
of colour basis and the strategy for the determination of the colour flows necessary for 
the coherent shower evolution are described in [28]. To the best of our knowledge, a large 
fraction of the matrix element calculations documented in this paper have never been 
performed before using other calculational tools. An independent implementation of the 
ALPHA algorithm has been exploited for the evaluation of multijet processes in hadronic 
collisions in [37, 6]. 

2.1.4 PDF sets and a s 

The code library includes a choice among some of the most recent PDF parameterizations. 
They can be selected at the beginning of the run through the variable ndns, which is 
mapped as follows 1 : 



ndns 


PDF 


[a s (m z )]n loop 


ndns 


PDF 


[a s (m z )]n loop 


1 


CTEQ4M [38] 


[0.116] 2 


101 


MRST99-1 [41] 


[0.1175] 2 


2 


CTEQ4L [38] 


[0.116] 2 


102 


MRST01-1 [42] 


[0.119] 2 


3 


CTEQ4HJ [38] 


[0.116] 2 


103 


MRST01-2 [42] 


[0.117] 2 


4 


CTEQ5M [39] 


[0.118] 2 


104 


MRST01-3 [42] 


[0.121] 2 


5 


CTEQ5L [39] 


[0.127]! 


105 


MRST01J [42] 


[0.121] 2 


6 


CTEQ5HJ [39] 


[0.118] 2 


106 


MRSTLO [43] 


[0.130]i 


7 


CTEQ6M [40] 


[0.118] 2 








8 


CTEQ6L [40] 


[0.118] 2 









In the case of NLO sets we use the 2-loop expression for a s 



a s (Q) = 



b 5 log(QVA|) (6 5 log 2 (QVA|))loglog(gVA|) 



(2.3) 



valid for Q > = 4.5 GeV, where 65 and 65 are the 1- and 2-loop coefficients of the 
QCD function, respectively, for 5 flavours. Threshold matching is applied in the case of 
Q < nib- In the case of LO sets (such as CTEQ*L or MRSTLO) we follow the prescriptions 



1 In addition to these default sets, we include in the package the full group of 40 CTEQ6M sets which 
allow the determination of PDF systematics errors [40] . To access these, the file alplib/ alppdf . f should be 
replaced in the linking stage with the file alplib/alppdf _cteq. f . A more complete set of parton densities, 
including old sets back to EHLQ and DO, is also available upon request. 
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used by the authors in performing the PDF fits. These vary from set to set. For exam- 
ple, set CTEQ5L was fitted using a LO expression for a s , while CTEQ6L used the NLO 
evolution. The table above lists the values of a s {mz) corresponding to the various sets, 
and indicates whether these values (and the relative evolution to different renormalization 
scales) correspond to the 1 or 2 loop formulation. 

2.1.5 Electroweak couplings 

In the current version of the ALPHA code the input couplings are derived from the standard 
SU (3) x SU{2) lxU (l)y tree level Lagrangian. The choice of input EW parameters deserves 
a short discussion. In ALPHA, the couplings of the electroweak (EW) and Higgs (H) bosons 
(including the respective selfcouplings) are parametrised in terms of the SU (2) coupling 
strength g, of the weak mixing angle sin of the electromagnetic fine structure constant 
a em , and of the masses of W, Z and H. We therefore have a total of 6 parameters needed 
to specify the value of general EW matrix elements. If we want to preserve gauge invariance 
we are however allowed to use only four independent parameters (plus fermion masses in 
the Yukawa sector). Treating mg as a free parameter, gauge invariance at the tree level 
demands that the remaining 5 parameters satisfy the following tree-level relationships: 

cos9 w = ^ (2.4) 
m z 

e = g sinOw . (2.5) 
The Higgs self-couplings and Yukawa couplings to fermions are furthermore given by: 

\ _9m 2 H 

Xhhh ~ w (2 - 6) 

2 2 

g L mjj 

Xhhhk = 32Ml <2 - ?) 

vr = -P- • (2-8) 

As a consequence, we cannot assign to the input parameters the values which are known 
from the current accurate experimental measurements, since these values are only consistent 
with the radiatively corrected versions of the above relations. This is not a merely formal 
issue: any tiny violation of the tree-level gauge relationships among the model parameters 
leads to violations of the equivalence theorem and leads to unphysical corrections to the 
tree-level cross-sections scaling like (Ey/My) 2n , n being the number of on/off-shell heavy 
gauge bosons appearing in the relevant diagrams and Ey, My their energy and mass 
respectively. In view of the large center of mass energy available in current and future 
hadron collisions, these spurious corrections could be numerically large, leading to the 
wrong high-energy behaviour of the cross sections. 

The current version of the code provides four choices for the setting of EW parameters. 
These choices are controlled by the variable iewopt, set at running time (default values for 
this variable are provided in the code, and are listed in the following sections describing the 
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individual hard processes). The different options are listed here below; the numerical values 
of the calculated parameters are obtained using the following set of inputs: mw = 80.41, 
m z = 91.188, sm 2 6 w = 0.231, a em (m z ) = 1/128.89, G F = 1.16639 x 10~ 5 : 



g = y/4:ira em (m z )/ sin 9 W = 0.6497, m w = g/\f 4\p2G F = 79.98, (2.9) 
mz = mw I cos 9w = 91.20 

iewopt=l. Inputs: mw, G F , sin 2 9w- Extracted: 

m z = mwl cos 6 w = 91.695, g = (4V2G F ) l/2 m w = 0.6532, (2.10) 
a e m(m z ) = {g sm9 w ) 2 /4ir = 1/127.51 

iewopt=2. Inputs: mz, a em {mz), sm 2 9w- Extracted: 

m w = m z cos 6 W = 79.97, g = \f^a em {m z )/ sm6 w = 0.6497 (2.11) 

iewopt=3. Inputs: mz, mw, Gf- Extracted: 

sin 2 W = 1 - {mw/mzf = 0.2224, g = (4V2G F ) 1/2 m w = 0.653, (2.12) 
aem(m z ) = (g sm8 w ) 2 /4ir = 1/132.42 . 

As a default, in all processes we employ iewopt=3. We verified that alternative options 
generate changes in the rates by at most few percent. Gauge and Higgs boson widths are 
calculated at tree level after the couplings have been selected. With the exception of the 
class of processes involving several gauge bosons, which will be discussed in detail Sec- 
tion 3.6 and where we set boson widths to 0, ALPHA uses fixed widths in the propagators. 

2.2 Unweighting and Shower evolution 

The starting point for the processing of events through the shower evolution is the gener- 
ation of a sample of unweighted events. This generation takes place through a two-step 
procedure (more details are given in the Appendix): 

1. to start with, a run of the parton-level code is performed as described in Section 2.1. 
Selecting the running mode option imode=l, weighted events are stored in a file. To 
limit the size of the file, instead of saving all the event information, we simply store 
the seed of the first random number used in the generation of the event, in addition 
to the event weight. 

2. at the end of the generation of the weighted event sample, the unweighting is per- 
formed by running the code once more, using a running mode option imode=2. In this 
running mode, the code will sequentially read the events stored in the file, and will 
perform the unweighting using the knowledge of the maximum weight of the sample, 



iewopt 



0. Inputs: a em (mz), G F , sin 2 9w- Extracted: 
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and of the weight of each individual event. When an event is selected by the unweight- 
ing, the seed of the random number is uploaded and all the information about the 
event (kinematics, flavours, spins and colours) is automatically reconstructed. The 
colour flow for the event is then selected according to the algorithm described in [28]: 
the subamplitudes corresponding to all colour flows compatible with the colour state 
of the event are first evaluated using ALPHA; one of them is then randomly extracted 
with a probability proportional to the squared modulus of the relative subamplitude. 
The momenta of the particles, together with their flavour and with the colour flow 
information, are then written to a file, which at the end of the run will contain the 
complete sample of unweighted events. 

At this point we are ready to process the events through the shower evolution. The stored 
events can be read by the chosen shower MC, the kinematics, colour and flavour information 
for each event being translated into the event format established by the Les Houches 
convention [34]. 

3. The available hard processes 
3.1 WQQ+ jets 

The physics content of the WQQ+ jets code has been described in detail in [28], where 
some phenomenological applications are also presented. We use the notation W as a short 
hand; what is actually calculated is the matrix element for a fermion-antifermion final 
state. All spin correlations and finite width effects are therefore accounted for. The quoted 
cross sections refer to a single lepton family; in the flavour assignment, the code selects 
by default an electron. Different flavours can be selected during the unweighting phase, 
covering all possible leptonic decays, as well as inclusive quark decays (for more details see 
the Appendix B.4). In the case Q = t, the top quark is left undecayed. The EW parameters 
are fixed by default using the option iewopt=3 (see eq. (2.10)). Only the leading-order EW 
diagrams are included in the calculation. 

The subprocesses considered include all configurations with up to 2 light-quark pairs; 
they are listed in Table 1, following the notation employed in the code. For all processes, 
the charge-conjugate ones are always understood. The above list fully covers all the possible 
processes with up to 3 light jets in addition to the heavy quarks. In the case of 4 extra 
jets, we do not calculate processes with 3 light-quark pairs. Within the uncertainties of 
the LO approximation, these can be safely neglected [4]. 

As a default, the code generates kinematical configurations defined by cuts applied to 
the following variables (the cuts related to the heavy quarks are only applied in the case 
of b, while top quarks are always generated without cuts): 

. PT v\ rf c \ AR 3j 

• PT b , r,\ AR bl 

• p T e , rf, p T u , AR ej , 
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jproc 


subprocess 


jproc 


subprocess 


jproc 


subprocess 


1 


qq' -► VFQQ 


2 


95 -► g'VFQQ 


3 


<?g -► g'VFQQ 


4 


ff5 qq'WQQ 


5 


gg' -► WQQq"q" 


6 


gg" -► WQQq'q" 


7 


g"g -► WQQq'q" 


8 


qq -► VFQQg'g" 


9 


qq' -► WQQgg 


10 


q'q -► VFQQgg 


11 


-> lUQQgg' 


12 


gg -► VFQQg'g 


13 


gg lUQQgg' 


14 


gg' WQQgg 


15 


gg' VFQQg'g' 


16 


g 5 - WQQq'q"q" 


17 


<?g - WQQq'q" q" 


18 


95 WQQqqq' 


19 


qg -> WQQq'qq 


20 


5g -► WQQqqq' 


21 


gq -► WQQq'qq 


22 


qg - WQQM? 


23 


gq - VFQQg'g'g' 


24 


gg - IUQQgg'g"g' 


25 


gg -> WQQqqqq' 











Table 1: Subprocesses included in the PFQQ+jets code. Additional final-state 
gluons arc not explicitly shown here but are included in the code if the requested 
light-jet multiplicity (AT < 4) exceeds the number of indicated final-state partons. 
For example, the subprocess jproc=l in the case of 2 light jets will correspond to 
the final state qq' — > WQQgg. The details can be found in the subroutine self lav 
of the file wqqlib/wqq. f . 



where AR a [> = y/[(r) a — T7&) 2 + (</> a — <f>b) 2 ]- The cut values can be provided by the user at 
run time. Additional cuts can be supplied by the user in the routine usrcut contained in 
the user file wqqwork/ wqqusr . f . 

In the code initialization phase, the user can select among 3 continuous choices for the 
parametrization of the factorization and renormalization scale Q: a real input parameter 
(qf ac) allows to vary the overall scale of Q, Q = qf ac x Qq, while the preferred functional 
form for Qq is selected through the integer input parameter iqopt: 



iqopt 





1 


2 


Ql 


m W + PTw 


m w 


m W + E m T 



where tut is the transverse mass defined as m\ = m 2 +p\, and the sum m\ extends 
to all final state partons (including the heavy quarks, excluding the W decay products). 

Some numerical benchmark results are given in Table 2 and 3. In the case of the 
minimal jet multiplicity, the results agree with previous calculations (see e.g. [44, 45]). 
The following scale and cuts are used: 

Q 2 = m 2 w +p 2 TW , (3.1) 

pr jct > 20 GeV, \r)j\ < 2.5, ARjj > 0.7 (3.2) 

p T b > 20 GeV, \rj b \ < 2.5, A% > 0.7, AR bj > 0.7 . (3.3) 

Here and in the following Sections we shall use the PDF set CTEQ5L. The quoted errors 
reflect the statistical accuracy of the integrations. We never tried to go beyond the percent 
level, to concentrate the CPU resources on the most computationally demanding channels. 
Results for the FNAL Tevatron refer to pp collisions at \fS = 2 TeV, those for the LHC 
refer to pp collisions at y/S = 14 TeV. 
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N = 


N = 1 


N = 2 


N = 3 


N = 4 


LHC (pb) 


2.222(4) 


3.013(9) 


1.83(1) 


0.831(8) 


0.307(5) 


FNAL (fb) 


332.2(7) 


86.2(4) 


18.3(2) 


3.17(3) 


0.44(3) 



Table 2: a(bb£v + N jets) at the Tevatron and at the LHC. Parameters and cuts 
are given in eqs. (3.1-3.3). 



LHC (fb) 


61.1(4) 


FNAL (fb) 


1.55(1) 



Table 3: u(tttv) at the Tevatron and at the LHC. 

3.2 Z/ 7 *QQ+ jets 

We use the notation Z/'y* as a short hand; what is actually calculated is the matrix element 
for a lepton-pair final state. All spin correlations and finite width effects are therefore 
accounted for. When the final state £ + £~ is selected, the interference between intermediate 
Z and 7* is also included. The quoted cross sections refer to a single lepton family; in 
the flavour assignement, the code selects by default the e + e~ pair. In the case of the final 
state vv the quoted cross sections include the decays to all 3 neutrino flavours, although 
we always label the neutrinos as v e . In the case Q = t, the top quark is left undecayed. 
The EW parameters are fixed by default using the option iewopt=3 (see eq. (2.11)). 

The subprocesses considered include all configurations with up to 2 light-quark pairs; 
they are listed in Table 4, following the notation employed in the code. For each process, the 
charge-conjugate ones are always understood. The above list fully covers all the possible 
processes with up to 3 light jets in addition to the heavy quarks. In the case of 4 or 
more extra jets, we do not calculate processes with 3 light-quark pairs. As in the case 
of associated production with a W, we expect that, within the uncertainties of the LO 
approximation, these can be safely neglected. 

As a default, the code generates kinematical configurations defined by cuts applied to 
the following variables (the cuts related to the heavy quarks are only applied in the case 
of b, while top quarks are always generated without cuts): 

• p T i et , rf ct , AR j:j 

• PT Q , V Q , ^Rqq 

• Pt £ , rf, ARfj, m(£ + £~), for £ + t~ final states 

• missing Et for vv final states. 

The cut on the dilepton invariant mass allows to optimise the sampling of the DY mass 
spectrum if the user is interested in events off the Z peak. Additional cuts can be supplied 
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jproc 


subprocess 


jproc 


subprocess 


jproc 


subprocess 


i 
i 


UU — > Z ljlj 


Z 


UU — > ZLJLJ 


q 
o 


gg — > Zt^t^ 


A 

4 


gu — > UZl^LJ 





„A , A7ru^ 

gu — > uzji^y^i 





Tic/ — > uZQQ 


7 


An , A7(~^C\ 

ug — > uzi^lj 




o 


gg — ► uuZQQ 


Q 

y 


„„ , AA7fS?^\ 

gg — > uuz: VcjVcj 


1 n 

1U 


uu — > uuZQQ 


1 1 
11 


A A s. AA7DH 

UU — ► UUZ Itfty 


1 9 
1Z 


uu — > uuZQQ 


i ^ 
io 


A A . AA7(lTS 

UU — ' UUZk^k^ 


1 /I 
14 


UU — > U U ZK^K^ 


10 


A A s. A 1 A 1 7DH 

UU — ► U U ZVcJi^ 


1 R 


UU — > UU Z/(^/(^/ 


1 7 
1 f 


UU — ► UU ZjVcJi^ 


15 


A A 1 v A A 1 7DH 

UU — ► UU Zv^Vej 


1 Q 


A A* , AAl 7(\~f} 

UU — > CLU Z/C^c^ 


9n 
zu 


„,„-; , AA7CU ; S 

UU — > UUZK^K^l 


91 
Zl 


AA . mm 7Cif} 

UU — > UUZK^v^l 


99 
ZZ 


„,A . r,.A7C\7) 

UU — ► UUZjk^K^I 


9Q 
Zo 


An, . Am 7 CU^ 

UU — > UUZLJLJ 


9/I 
Z4 


„,A s. mA 7 fl~f\ 

UU — > UUZv^l^) 


9^ 
ZD 


Am . Am 7fS7} 

UU — > UUZi^i^ 


9fi 
ZO 


UU — ► 00ZLJLJ 


97 
Z / 


A A i. hh 7DH 

UU — > OOZLjLJ 


zo 


gg — ► oo Zj LJLJ 


on 

zy 


gu — > uuuZQQ 


OA 

oil 


ug — ► uuuZQQ 


31 


gu — ► uu'u'ZQQ 


32 


ug — ► uu'u'ZQQ 


33 


gu — > uddZQQ 


34 


ug — ► uddZQQ 


35 


gu — > ubbZQQ 


36 


ug — > ubbZQQ 


37 


gd — > dddZQQ 


38 


dg — > dddZQQ 


39 


gd -► dd'd'ZQQ 


40 


dg dd'd'ZQQ 


41 


gd — > duuZQQ 


42 


dg — ► duuZQQ 


43 


gd -► dbbZQQ 


44 


dg dbbZQQ 







Table 4: Subprocesses included in the Z/7*QQ+jets code. It is always understood 
that quarks u and c? represent generic light quarks of type up or down. The Z in 
the table stands for a neutral £ + £~ (yv) lepton pair. Additional final-state gluons 
are not explicitly shown here but are included in the code if the requested light- 
jet multiplicity (N < 4) exceeds the number of indicated final-state partons. For 
example, the subprocess jproc=l in the case of 2 light jets will correspond to the 
final state uu — > ZQQgg. The details can be found in the subroutine self lav of the 
file zqqlib/zqq. f . 

by the user in an appropriate routine. The choice of factorization and renormalization scale 
is similar to what given for the WQQ processes, with the W mass replaced by the mass 
of the DY pair (if the DY mass range excludes the value of mz), or by mz (when the DY 
mass range includes mz)- 

Some benchmark results are given in Table 5, obtained for the following set of cuts 
and scale choice: 

Q 2 = m \ + p\ z , 80 GeV < m u < 100 GeV (3.4) 
p r jct >20GeV, \m\<2.5, AR j3 > 0.7 (3.5) 
p T b > 20 GeV, \r] b \ < 2.5, A% > 0.7, R bj > 0.7 . (3.6) 

3.3 W+ jets 

As in the previous cases, we use the notation W as a short hand; what is actually calculated 
is the matrix element for a lepton+neutrino final state. All spin correlations and finite 
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N = 


TV = 1 


N = 2 


N = 3 


N = 4 


LHC, (fb) 


1840(5) 


1085(4) 


444(3) 


154(2) 


44(1) 


FNAL, (fb) 


49.3(1) 


13.18(5) 


2.57(2) 


0.400(4) 


0.0511(5) 



Table 5: <r(£ + £ bb + N jets) at the Tevatron and at the LHC. Parameters and cuts 
are given in eqs. (3.4-3.6). 



jproc 


subprocess 


jproc 


subprocess 


jproc 


subprocess 


1 


gg' -»• W 


2 


qg -> g'PF 


3 


gg -> g'PF 


4 


gg -> qq'W 


5 


gg' lFg"g" 


6 


gg" -► JFg'g" 


7 


q"q -► JFg'g" 


8 


gg — > VFg'g" 


9 


gg' — > VFgg 


10 


g'g — > IFgg 


11 


gg — > Wqq 1 


12 


gg — > IFg'g 


13 


— > VFgg' 


14 


qq' — > VFgg 


15 


gg' — > TFg'g' 


16 


qg -► IFg'g"g" 


17 


gg — > Wq'q"q" 


18 


gg -> IFggg' 


19 


95 VFg'gg 


20 


gq -> IFggg' 


21 


gg -> IFg'gg 


22 


g5 -> Wq'q'q' 


23 


gg -> Wq'q'q' 


24 


gg -> Wqq'q"q' 


25 


55 -> VFgggg' 











Table 6: Subprocesses included in the VF+jets code. Additional final-state gluons 
are not explicitly shown here but are included in the code if the requested light- 
jet multiplicity (N < 6) exceeds the number of indicated final-state partons. For 
example, the subprocess jproc=l in the case of 2 jet will correspond to the final 
state qq' — > Wgg. The details can be found in the subroutine self lav of the file 
wjetlib/wjets . f . 

width effects are therefore accounted for. The quoted cross sections refer to a single lepton 
family; in the flavour assignment, the code selects by default an electron. Different flavours 
can be selected during the unweighting phase, covering all possible leptonic decays, as well 
as inclusive quark decays (for more details see the Appendix B.4). The EW parameters 
are fixed by default using the option iewopt=3 (see eq. (2.10)). 

The subprocesses considered include all configurations with up to 2 light-quark pairs; 
they are listed in Table 6, following the notation employed in the code. For each process, 
the charge-conjugate ones are always understood. The above list fully covers all the possible 
processes with up to 3 light jets in addition to the heavy quarks. In the case of 4 extra 
jets, we do not calculate processes with 3 light-quark pairs. Within the uncertainties of 
the LO approximation, these can be safely neglected [4]. 

As a default, the code generates kinematical configurations defined by cuts applied to 
the following variables: 

. p T Q, r,* AR Q q 
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N = 


N = 1 


N = 2 


N = 3 


N = 4 


N = 5 


iV = 6 


LHC (pb) 


18068(4) 


3412(4) 


1130(2) 


342.9(1.4) 


100.6(1.4) 


27.6(4) 


7.14(15) 


FNAL (pb) 


2087.0(6) 


225.8(2) 


37.3(2) 


5.66(6) 


0.745(4) 


0.0864(15) 


0.0086(2) 



Table 7: er(^ + TV jets) at the Tevatron and at the LHC. Parameters and cuts are 
given in cqs. (3.7-3.8). 



• p T , V , PT U , ARij ■ 

The respective threshold values can be provided by the user at run time. Additional cuts 
can be supplied by the user in an appropriate routine. The choice of scale follows the same 
conventions as for the WQQ case. 

Some benchmark results are given in Table 7, obtained for the following set of cuts 
and scale choice: 

Q 2 = m 2 w +p 2 TW , (3.7) 
Pr jct >20GeV, \r}j\ < 2.5, ARjj > 0.7 . (3.8) 

For processes with up to 4 jets, we verified the numerical agreement with the results 
of the VECBOS code [4]. 

3.4 W + c+ jets 

The processes included in this code are a subset of the ones treated in section 3.3. Here 
the final states consisting exclusively of one c (or c) quark in association with a W and 
additional light jets are singled out. Events with charm quark pairs are not included, and 
should be generated using the wqq processes. All spin correlations and finite width effects 
in the fermionic decay of the W are accounted for. The W decay modes can be selected 
when running the code with imode=2 to produce unweighted events, as discussed in the 
appendix B.4. The EW parameters are fixed by default using the option iewopt=3 (see 
eq. (2.10)). 

The subprocesses considered include all configurations with up to 2 light-quark pairs 
(where "light" includes the charm); they are listed in Table 8, following the notation 
employed in the code. For each process, the charge-conjugate ones are always understood. 
The above list fully covers all the possible processes with up to 2 light jets in addition to the 
heavy quarks. In the case of 3 extra jets, we do not calculate processes with 3 light-quark 
pairs. As in the case of VF+jet production, we expect these to be negligible. 

As a default, the code generates kinematical configurations defined by cuts applied to 
the following variables: 

• PT i ct , if ct , ARjj 

• PT e , V e , PT U , ARij , 
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jproc 


subprocess 


jproc 


subprocess 


1 


/ XT T 

gd -> Wc 


2 






3 


gg — > IFcc' 


4 


gg" 


-► Wq'q"{q'/q" = c 


5 


g"g Wq'q"(q'/q" = c) 


6 


-► Wq'q"(q'/q" = c) 


7 


qq -> Wqq'(q/q' = c) 


8 


gg - 


-> Wq'q(q/q' = c) 


9 


qq — ► VFgg'(g/g' = c) 


10 




-+ VFcg"g" 


11 


5 c' -► WcgV 


12 




-► VFc'c'c 


13 


c'g -► VFcc'c' 


14 




-+ VFc'c'c 


15 


gc' — > Wcc'd 


16 






17 


gg — ► Wd(tdc 









Table 8: Subprocesses included in the W + c+jets code. The symbol d refers to 
either of the two Cabibbo partners of the charm quark, d or s. Additional final-state 
gluons arc not explicitly shown here but are included in the code if the requested 
light-jet multiplicity (N < 5) exceeds the number of indicated final-state partons. 
For example, the subprocess jproc=l in the case of 2 jet will correspond to the 
process gd — > W egg. The details can be found in the subroutine self lav of the file 
wcjetlib/wcjets .f . 





N = 


N = 1 


N = 2 


N = 3 


N = 4 


N = 5 


LHC (pb) 


418.9(8) 


183.7(8) 


55.6(3) 


14.9(1) 


3.65(3) 


0.848(7) 


FNAL (fb) 


8740(20) 


2390(10) 


360(2) 


42.2(2) 


4.14(2) 


0.353(2) 



Table 9: o{lv + c + N jets) at the Tevatron and at the LHC. Parameters and cuts 
are given in eqs. (3.9-3.10). 



where the cuts on c-quarks are the same as for light jets. The respective threshold values 
can be provided by the user at run time. Additional cuts can be supplied by the user in 
an appropriate routine. The choice of scale follows the same conventions as for the WQQ 
case. Some benchmark results are given in Table 9, obtained for the following set of cuts 
and scale choice: 

Q 2 = m 2 w +p 2 TtW , (3.9) 
p T ict > 20 GeV, \r)j\ < 2.5, ARjj > 0.7 . (3.10) 

3.5 Z/ 7 *+ jets 

We use the notation Z/j* as a short hand; what is actually calculated is the matrix element 
for a charged lepton or neutrino pair final state. All spin correlations and finite width effects 
are therefore accounted for. When the final state £ + £~ is selected, the interference between 
intermediate Z and 7* is also included. The quoted cross sections refer to a single lepton 
family; in the flavour assignment, the code selects by default the e + e~ pair. In the case of 
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the final state vv the quoted cross sections include the decays to all 3 neutrino flavours, 
although we always label the neutrinos as v e . The EW parameters are fixed by default 
using the option iewopt=3 (see eq. (2.11)). 

All subprocesses with up to 2 light quark pairs are included. This means that the 
cross-sections with up to 3 final-state partons are exact. The emission of additional hard 
gluons can however be calculated, and the current version of the code works with up to a 
6 final-state jets. The subprocesses considered are listed in Table 10. 



jproc subprocess 



1 


uu 


-► Z 


4 


gd 


-► dZ 


7 


99 


— > uuZ 


10 


dd 


— > ddZ 


13 


uu 


— > u'u'Z 


16 


uu' — > uu' Z 


19 


uu 


-► ddZ 


22 


du 


— > udZ 


25 


gu 


— > uuuZ 


28 


ug 


— ► uu'u'Z 


31 


gd 


— > dddZ 


34 


dg 





jproc 


subprocess 


2 


Z 


5 


u<7 — > 


8 


<7<7 — ► ddZ 


11 


mu — ► miZ 


14 


-► d'd'Z 


17 


dd' -► dd'Z 


20 


dd — ► miZ 


23 


ud — ► udZ 


26 


ug — > uuuZ 


29 


gu — > uddZ 


32 


dg — ► dddZ 


35 


gd — > duuZ 



jproc 


subprocess 


3 


git — ► uZ 


6 


d 5 -► dZ 


9 


tm — ► traZ 


12 


dd — > ddZ 


15 


uu' — > u?/Z 


18 


dd' -► dd'Z 


21 


ud — ► itdZ 


24 


du — > duZ 


27 


gu — > uu'u'Z 


30 


ug — > uddZ 


33 


gd dd'd'Z 


36 


d</ — > duuZ 



Table 10: Subprocesses included in the Z/7*+jets code. It is always understood 
that quarks u and d represent generic quarks of type up or down. The Z in the table 
stands for a neutral £ + £~ (vv) lepton pair. The complex conjugate processes are 
also understood. Additional final-state gluons are not explicitly shown here but are 
included in the code if the requested jet multiplicity (N < 6) exceeds the number of 
indicated final-state partons. For example, the subprocess jproc=l in the case of 2 
jets will correspond to the final state uu — > Zgg. The details can be found in the 
subroutine self lav of the file zjetlib/zjets . f . 



As a default, the code generates kinematical configurations defined by cuts applied to 
the following variables: 

• P T ict , rf ct , ARjj 

• pt £ , t) , ARij for final states 

• missing Et for vv final states. 

Additional cuts can be supplied by the user in an appropriate routine. The choice of scale 
follows the same conventions as for the ZQQ case. 

Some benchmark results are given in Table 11, obtained for the following set of cuts 
and scale choice: 

Q 2 = m 2 z +p 2 TZ , (3.11) 
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N = 


N = 1 


N = 2 


N = 3 


N = 4 


N = 5 


N = 6 


LHC (pb) 


1526(1) 


320.9(5) 


104.6(2) 


31.6(2) 


9.4(2) 


2.51(4) 


0.65(2) 


FNAL (pb) 


179.4(2) 


21.44(2) 


3.36(1) 


0.489(2) 


0.0630(3) 


0.00700(4) 


0.000690(6) 



Table 11: a(£ + £ + N jets) at the Tevatron and at LHC. Parameters and cuts are 
given in cqs. (3.7-3.8). 



p T ict > 20 GeV, 1^1 < 2.5, AR j:j > 0.7 (3.12) 
80 GeV < mu < 100 GeV . (3.13) 

3.6 nW + mZ + 1H + k jets 

The code computes all processes where massive EW gauge bosons and/or Higgs particles 
are produced on shell plus (up to 3) additional light jets. The limitations in the number of 
final state particles is therefore the following: k < 3 and n + m + l + k<8. It is intended 
that at least one EW gauge boson or Higgs particle appear in the final state (n+m + l > 0). 
All contributions from QCD and EW processes are included. In particular, all processes of 
the gauge-boson-fusion type are present when the number of final state quarks is at least 
2 (e.g. qq — ► qqH). All gauge bosons are decayed to fermion pairs, taking therefore into 
account all spin correlations among the decay products by means of exact matrix elements. 
More information on gauge-invariance issues in presence of decays are discussed below, and 
details on the choice of the decay modes can be found in Appendix B.4. The Higgs boson 
decays will be soon available. 

In the case of m = I = (n = I = 0) the code reproduces the results of the pro- 
grams described in Sections 3.3 and 3.5, up to the following effects: final-state finite-width 
corrections are absent here since the gauge bosons are kept on shell; branching ratios de- 
pend on the selected decay mode for the Z, while they are kept to 1 in the case of the W 
bosons, whose final states are selected during the unweighting phase (see Appendix B4); 
the 7* — ► £ + £~ contributions are only present in the code of Section 3.5; EW production of 
jets (W — ► jj, Z/j* — > jj, plus EW boson exchanges between quarks) are included here. 

The EW parameters are fixed by default using the option iewopt=3 (see eq. (2.12)). In 
this way we are guaranteed that all gauge boson masses are consistent with the measured 
values. The coupling strengths extracted from the fixed inputs, however, will differ by few 
percents from their best values. The user can select the option iewopt=0 (see eq. (2.9)), 
where the couplings are matched to the radiatively corrected best values, at the price of 
working with gauge boson masses which differ by few percents from the measured masses. 
Both schemes are equally consistent at the LO. We verified that cross-sections obtained 
using the two schemes differ from each other at the level of few percent, a negligible effect 
compared to the large uncertainties related to the choice of factorization and renormaliza- 
tion scales. 
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n odd 



jproc 


subprocess 


jproc 


subprocess 


jproc 


subprocess 


1 


Q a> _». jr± 


2 


00 — > Q 1 J~ 


3 


oo — > o ' J- 


4 


99 ^'-^ 


5 


qq' -► JF±9"g" 


6 


qq" — ► T^q'q' 


7 


g"9 -► .FVg" 


8 


g"g" -► ^±99' 


9 


qq' — ► T^qq 


10 


g'g -> JT±gg 


11 


09 ^99' 


12 


qq -> .FVg 


13 


gg -> ./^gg' 


14 


99' -> T±qq 


15 


99 — > 9 9 


16 




17 


gq - ^yg" 


18 


95 -> .^±999' 


19 


gg -> J r± q'qq 


20 


59 -> .T^ggg 7 


21 


59 -> T^cfqq 








n even 







jproc 


subprocess 


jproc 


subprocess 


jproc 


subprocess 


1 


99^° 


2 


59 9^° 


3 


95 -> 9^"° 


4 


55 -> 99^ 


5 


99 -> 99^° 


6 


99 — ► q'q'T 


7 


99' — > qq'J 70 


8 


99" — > qq" 


9 


qq' 9 ,/ 9 /,/ ^° 


10 


qq -> 99^° 


11 


qq — > q'q' J 70 


12 


99 — > 99 > 


13 


99' — ► 99' or qq'J 70 


14 


qq" -> gg".? 70 


15 


99' -» gV'-F 


16 


qq" g'g'"^ 


17 


59 -> 999^° 


18 


95 -> 999^° 


19 


59 -> qq'q'F 


20 


95 -> qq'q'F 


21 


59 - 9?V^ 


22 


qg - gg"g"^° 


23 


gq — > q'q'qJ 70 


24 


95 -> g'g'g.? 70 


25 


59 -► q'q"q'"F° 


26 


qg q'q"q'"T° 







Table 12: Subprocesses included in the nW^ + mZ + liY+jets code, for the cases 
n=odd (nVK + mZ + 111 = ^) and n=cvcn (nW + mZ + IH = It is always 

understood that quarks g and g' belong to the same iso-doublct, while q and g" belong 
to different iso-doublets. Additional final-state gluons are not explicitly indicated but 
are included in the code. The details can be found in the subroutine self lav of the 
file vbjetslib/vbjets.f . 

The subprocesses considered are listed in Table 12, where a distinction is made ac- 
cording to the number of final state W's. If n is odd, the whole Cabibbo structure of the 
matrix element is correctly taken into account. If n is even, we work in the cos# c = 1 
approximation. The reason of this choice is that, in the latter case, different Cabibbo 
structures may interfere at the amplitude level. The quantity T^- in Table 12 stands for 
nW + mZ + IH with n odd, while J^ stands for nW + mZ + IH when n is even. 

As a default, the code generates kinematical configurations denned by cuts applied to 
the following variables: 

• p T i ct , rf ct , ARjj . 

Additional cuts can be supplied by the user in an appropriate routine. Cuts on the lep- 
tons or on the Higgs decay products should be imposed in the analysis routine, using the 
momentum and flavour information as specified in Appendix B.4. 
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In the code initialization phase, the user can select among 3 continuous choices for the 
parametrization of the factorization and renormalization scale Q: a real input parameter 
(qf ac) allows to vary the overall scale of Q, Q = qf ac x Qo, while the preferred functional 
form for Qo is selected through an integer input parameter (iqopt=0,l,2). In more detail: 



iqopt 





1 


2 


Qo 


(M B ) 







where (Mb) and Yl Mb are the average and sum of the EW boson masses (vectors as 
well as Higgses), respectively. 

Some benchmark results are given in Tables 13-18, using the scale Q = J^Mb, with 
the sum extended over all bosons in the final state. 

In order to deal, in a gauge invariant way, with the problem of the unstable bosons 
appearing in the intermediate states, we computed the matrix element by using zero-width 
propagators and by cutting away events around the mass of the unstable particle M in 
such a way to keep the area of a Breit-Wigner distribution. In other words, for a particle 
of mass M and width V, our effective cut so is determined by the equation 

M 2 ~so i rM 2 ^ 

„ & FW = L*Fm¥Tw' (3 - 14) 

giving so = 2 ™ ■ In practice the program discards events for which |s — M 2 \ < winsize x 
rM, with the adjustable parameter winsize set, by default, to the value winsize = ^. 
The described procedure gives sensible results when the portion of resonance cut away is of 
the order of a few GeV. Otherwise, holes start becoming visible in the distributions. This 
is safe for Z and W vector bosons. But problems arise for very heavy Higgses. For this 
reason we put a protection in the code to inhibit the calculation when Tjj > 10 GeV. 

We checked that our algorithm reproduces, within few percents, the results obtained 
by using the following approach: we set T = in the matrix element and multiply the final 
result by the factor 



As a consistency check of our calculations for large multiplicities of gauge bosons, 
we verified that the production rates for multiple gauge bosons when the Higgs mass is 
above the threshold for diboson decay are well approximated by the incoherent sum of 
the processes mediated by an on-shell Higgs, plus those where the Higgs contribution is 
suppressed in the intermediate states. In the tables we show, this is for example seen in 
the comparison of the WWjj vs Hjj or VVV vs HV rates at run =200 GeV. In the case 
of 3 gauge boson production, we verified the agreement with the results shown in ref. [46], 
up to the WZZ channel, whose rate is erroneously reported in the tables of Section 5.33 2 . 
In the case of four boson production, a previous calculation has been documented in [47]. 
A comparison of the results of our code with the numbers presented in that paper shows 
however some discrepancies. 

2 A. Ghinkulov, private communication. 
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WW 


WZ 


ZZ 


LHC (pb) 


74.8(5) 


28.1(2) 


10.85(6) 


FNAL (pb) 


8.51(5) 


2.45(1) 


1.027(5) 



Table 13: Diboson production at the Tevatron and at the LHC. 



WW+}ets 


jj, central 


jjj, central 


jj, fwd 


jjj, fwd 


LHC, m H =120 (pb) 


18.6(2) 


8.16(9) 


0.256(4) 


0.365(9) 


LHC, m H =200 (pb) 


18.9(2) 


8.3(1) 


0.546(24) 


0.389(9) 


FNAL, m H =120 (fb) 


336(1) 


49.1(2) 


0.201(1) 


0.0789(3) 


FNAL, m H =200 (fb) 


364(2) 


54.9(8) 


0.415(4) 


0.096(2) 



Table 14: Associated production of WW and jets, at the Tevatron and at the LHC. 
In all cases, Exj > 20 GeV and ARjj > 0.7. The central configurations correspond 
to all jets with \r]j\ < 2.5. The fwd configurations have two jets in opposite rapidity 
hemispheres with 2.5 < \rjj\ < 5, plus a central jet in the jjj case. 



3.7 QQ+ jets 

The subprocesses considered include all configurations with up to two light quark pairs. 
They are listed in Table 19, following the notations employed in the code. The list covers all 
possible processes involving up to 4 light-parton jets with 2 light quark pairs. Subprocesses 
involving 3 light-quark pairs, which would only appear in the case of 4 jets in addition to 
the heavy quarks, are not included, as they are expected to contribute a negligible rate. 

As a default, the code generates kinematical configurations defined by cuts applied to 
the following variables (the cuts related to the heavy quarks are only applied in the case 
of b, while top quarks are always generated without cuts): 

. PT ^, if et , ARjj^Rjt 

• PT b , 7] b , A% . 

The respective threshold values can be provided by the user at run time. Additional cuts 
can be supplied by the user in an appropriate routine. In the case Q = t the decay of 
ti pairs in six fermions (in the on-shell approximation) is enforced. With this option the 
decay of the top quarks takes into account all spin correlations among the decay products 
by means of exact matrix elements. More information on the selection of decay products 
can be found in Appendix B.4. 

In the code initialization phase, the user can select between 2 choices for the parame- 
terization of the factorization and renormalization scale Q. A real input parameter (qf ac) 
allows to vary the overall scale of Q, Q = qf ac x Qq, while the preferred functional form 
for Qq is selected through an integer input parameter (iqopt=0,l). In more detail: 
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H+}ets 


jj, central 


j j j, central 


■ " C 1 

jj, fwd 


• • • p 1 

JJJ,fwd 


LHC, mn =120 (pb) 


1.27(1) 


0.458(6) 


0.504(3) 


0.089(1) 


LHC, m# =140 (pb) 


0.96(1) 


0.320(4) 


0.458(2) 


r\ n'70 / 

0.078(2) 


LHC, m H =200 (pb) 


0.459(4) 


0.132(3) 


0.346(2) 


0.0528(6) 


FNAL, =120 (fb) 


121(1) 


26.7(2) 


0.769(2) 


0.0730(5) 


FNAL, m H =140 (fb) 


75.0(4) 


16.0(1) 


0.573(1) 


0.0513(3) 


FNAL, m H =200 (fb) 


21.8(2) 


4.28(3) 


0.246(1) 


0.0193(1) 



Table 15: Associated production of H and jets, at the Tevatron and at the LHC. 
In all cases, Exj > 20 GeV and ARjj > 0.7. The central configurations correspond 
to all jets with \t]j\ < 2.5. The fwd configurations have two jets in opposite rapidity 
hemispheres with 2.5 < \r]j\ < 5, plus a central jet in the jjj case. 





WWW 


WWZ 


WZZ 


zzz 


LHC, m H =120 


130(1) 


98(2) 


31.0(4) 


10.9(1) 


LHC, m H =200 


305(5) 


199(4) 


95(1) 


45.2(5) 


FNAL, m H =120 


6.39(2) 


5.13(4) 


1.21(1) 


0.519(2) 


FNAL, m H =200 


18.1(2) 


13.5(2) 


5.47(6) 


3.36(3) 



Table 16: Triboson production at the Tevatron and at the LHC (fb). 



iqopt 





1 


Qo[tt? 


o 

mf 


m\ + (p2 ) 


Qo[bb} 2 


s 


(Pi) 



where (pr) is the average pt of light partons, if Q = t, while is the average transverse 
momentum of all final state jets in the case of Q = b. 

Some benchmark results are given in Table 20, obtained for the following set of cuts 
and scale choice: 

it: m t = 175 GeV, Q 2 = m\ (3.15) 
bb: m b = 4.75 GeV, Q 2 = (p 2 Tb + P^, + J>!>>/(2 + AO (3-16) 
p T ict > 20 GeV, \T)j\ < 2.5, ARjj > 0.7 (3.17) 
p T b > 20 GeV, \r) b \ < 2.5, AR bl > 0.7, AR bj > 0.7 . (3.18) 

In the case of and 1 jet, we find agreement with the results obtained using the 0(o$.) 
code of ref . [48] . 
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WWWW 


WWWZ 


WWZZ 


WWWWW 


WWWWZ 


LHC, m H =120 


0.606(6) 


0.72(1) 


0.48(1) 


5.8(1) • 10" 3 


10.8(3) • 10~ 3 



Table 17: 4- and 5-boson production at the LHC (fb). 



WH 


m H =120 


m H =140 


m H =200 


ZH 


m H =120 


m H =140 


m H =200 


LHC (pb) 


1.364(8) 


0.833(4) 


0.251(2) 




0.727(4) 


0.449(2) 


0.137(1) 


FNAL (fb) 


122.4(6) 


70.0(2) 


16.55(4) 




75.0(3) 


44.2(2) 


11.24(3) 



Table 18: H+boson production at the Tevatron and LHC. 



jproc 


subprocess 


jproc 


subprocess 


jproc 


subprocess 


1 


99 QQ 


2 


gg -> QQ 


3 


3g -> QQg 


4 


qg -> QQg 


5 




6 


gg -> QQgg 


7 


-> QQg'g' 


8 


gg' -> QQgg' 


9 


gg' -»■ QQgg' 


10 


gg -> QQgg 


11 


gg -> QQqq'q' 


12 


qg -> QQqq'q' 


13 


gg -> QQggg 


14 


qg -> QQggg 


15 


99 QQqqqq 


16 


<7<7 ->• QQqqq'q' 











Table 19: Subprocesses included in the QQ+jets code. Additional final-state gluons 
are not explicitly shown here but are included in the code if the requested light- 
jet multiplicity (N < 6) exceeds the number of indicated final-state partons. For 
example, the subprocess jproc=l, in the case of 2 extra jets, will correspond to the 
final state gg — > QQgg. The details can be found in the subroutine self lav of the file 
2Qlib/2Q.f. For each process, the charge-conjugates ones are always understood. 

3.8 QQQ'Q'+ jets 

The list of available processes is given in Table 21. This covers all possible processes up to 
2 light jets and at most 1 light quark pair. The cases with one additional heavy quark pair 
(e.g. tibbbb) are also included. Subprocesses with two light-quark pairs give a negligible 
contribution to the case of 2 extra jets, and are not calculated. For the parameterization 
of the factorization and renormalization scale the user has two choices, iqopt=0,l, as 
described in the following Table: 



iqopt 





1 


Qo[bb] 2 


s 


(P 2 T ) 


Q [bt] 2 


s 


2 


Q [tt] 2 


s 


mi 



where {p?n) is the average p 2 -. of light partons and 6's. 
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QQ + N jets 


N = 


N = 1 


iV = 2 


N = 3 


N = 4 


N = 5 


N = 6 


Q = t, LHC (pb) 


530.0(8) 


462.6(6) 


255(1) 


111.5(5) 


42.4(4) 


14.07(16) 


4.36(8) 


Q = t, FNAL (fb) 


6,364(8) 


1,592(3) 


282(1) 


40.6(3) 


4.83(4) 


0.483(6) 


0.0419(9) 


Q = b, LHC (nb) 


1,533(4) 


422(1) 


130.2(6) 


30.9(4) 


7.5(4) 


1.53(6) 


0.337(9) 


Q = b, FNAL (pb) 


72,1(1) 


12,15(2) 


2,51(1) 


365(7) 


47.3(9) 


5.6(2) 


0.58(2) 



Table 20: a(QQ + N jets) at the Tevatron and at the LHC, with parameters and 
cuts given in eqs. (3.15-3.18). 



jproc subprocess 
1 gg - QQQ'C? 

4 g<7 QQQ'Q'q 



jproc subprocess 
2 QQQ'Q' 

5 55 ->• QQQ'Q qq 



jproc subprocess 

3 gq -► QQQ'Q'q 

6 5 5 QQQ'Q'bb 



Table 21: Subprocesses included in the QQQ'Q'+ jets code. Additional final-state 
gluons are not explicitly shown here but are included in the code if the requested 
light-jet multiplicity (N < 4) exceeds the number of indicated final-state partons. 
For example, the subprocess jproc=l, in the case of 2 extra jets, will correspond 
to the final state gg — ► QQQ'Q gg. The details can be found in the subroutine 
self lav of the file 4Qlib/4Q.f. For each process, the charge-conjugates ones are 
always understood. 



Again, as a default, the code generates kinematical configurations defined by cuts 
applied to the following variables (the cuts related to the heavy quarks are only applied in 
the case of b, while top quarks are always generated without cuts): 

. p T J ct , ??i ct , ARjj,AR jb 

The respective threshold values can be provided by the user at run time. Additional cuts 
can be supplied by the user in an appropriate routine. 

Some benchmark results are given in Table 22, obtained for the following set of cuts 
and scale choice: mt=175 GeV, mfe=4.75 GeV, 

tttt and tibb : Q 2 = m\ (3.19) 

bbbb : Q 2 = (p 2 Ttb +p 2 Ttl + J>^)/(2 + N) (3.20) 

p r jct >20GeV, \rjj\ < 2.5, AR j3 > 0.7 (3.21) 

p T b > 20 GeV, \rn,\ < 2.5, AR bi > 0.7, AR bj > 0.7 . (3.22) 
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QQQQ + iV jets 


iV = 


N = 1 


N = 2 


N = 3 


N = 4 


iiit, LHC (id) 


12.73(8) 


17.4(2) 


13.5(1) 


7.55(6) 


3.48(5) 


tt"66, LHC (pb) 


1.35(1) 


1.47(2) 


0.94(2) 


0.457(8) 


0.189(4) 


tt"66, FNAL (fb) 


3.44(3) 


0.95(1) 


0.154(1) 


0.0187(2) 


0.00187(5) 


bbbb, LHC (pb) 


477(2) 


259(5) 


95(1) 


28.6(6) 


25.0(3) 


bbbb, FNAL (pb) 


6.64(5) 


2.25(3) 


0.470(5) 


0.076(1) 


0.0025(5) 



Table 22: a(QQQ'Q' + N jets) at the Tevatron and at the LHC, with parameters 
and cuts given in eqs. (3.19-3.22). 



jproc 


subprocess 


jproc 


subprocess 


jproc 


subprocess 


1 


gg - QQtf 


2 


gg QQif 


3 


#g -> (3Qg# 


4 


95 QQqH 


5 


5ff QQqqH 


6 




7 


qq — ► QQqqH 


8 


qq' QQgg'tf 


9 


— > QQqqH 


10 


qq' -► QQgg'if 


11 


QQq'q'H 


12 


qq -► QQtotf 


13 


c/g -> QQqqqH 


14 


qg -> QQqqqH 


15 


#g -> QQq'q'qH 


16 


qg -> QQq'q'qH 


17 


gq — > QQbbqH 


18 


qg — > QQbbqH 



Table 23: Subprocesses included in the QQH code. The details can be found in the 
subroutine self lav of the file QQhlib/QQh . f . For each process, the charge-conjugate 
subprocesses are always understood. 



3.9 QQH+ jets 

The list of processes is given in Table 23. The Higgs is produced only via Yukawa couplings 
to the heavy quarks, no other EW process is included. All cases with up to 2 light-quark 
pairs are included, covering in full the possible final states with up to 3 jets in addition to 
the QQH system. The code will otherwise deal with up to 4 extra jets. 

As a default, the code generates kinematical configurations defined by cuts applied to 
the following variables: 

. PT ^\ rf et , ARjj,AR jb 

• PT b , r] b , AR bi . 

The respective threshold values can be provided by the user at run time. Additional cuts 
can be supplied by the user in an appropriate routine. In the case Q = t the decay of tt 
pairs in six fermions (in the on-shell approximation). With this option the decay of the 
top quarks takes into account all spin correlations among the decay products by means of 
exact matrix elements. More information on how to use the selection of decay products 
can be found in Appendix B.4. The Higgs boson decay will be soon available. 
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QQH 


rriH = 120 GeV 


rriH = 150 GeV 


m,H = 200 Gev 


Q = t, LHC (fb) 


401(2) 


212(1) 


89.1(4) 


Q = t, FNAL (fb) 


4.22(2) 


2.00(1) 


0.637(4) 


Q = b, LHC (fb) 


599(3) 


279(3) 


99(2) 


Q = b, FNAL (fb) 


3.73(3) 


1.20(1) 


0.240(2) 



Table 24: <r(QQH), at the Tevatron and at the LHC, with parameters given in 
eqs. (3.23-3.24). No cuts applied. 



For the parameterization of the factorization and renormalization scale the user has 
three choices, iqopt=0,l,2, as described in the following Table: 



iqopt 





1 




Qo[bb} 2 


s 


m\ + (pi) 


(P 2 T ) 


Qo[tt} 2 


s 


(2m t + m H ) 2 


(P 2 T ) 



where (p^) is the average px 2 of light and heavy partons. Some benchmark results are 
given in Table 24, obtained for the following set of cuts and scale choice: 

it: Q 2 = {2m t + m H f (3.23) 
bb: Q 2 = m 2 H + ( P 2 Tb+P 2 Trb )/2. (3.24) 

As a default, the Yukawa coupling of the Higgs to the bottom quarks is evaluated at the 
b-quark pole mass. Since the rate is directly proportional to y 2 , the result corresponding 
to the choice yi, oc m^mx) can be obtained via a trivial rescaling. The numbers we obtain 
in the case of extra jets agree with what found in the literature [49, 50], after possibly 
correcting for the difference between pole and running 6-mass. 

3.10 N jets 

The subprocesses considered include all configurations with up to two light quark pairs. 
They are listed in Table 25, following the notations employed in the code. The list covers all 
possible processes involving up to 6 light-parton jets with 2 light quark pairs. Subprocesses 
involving 3 light-quark pairs, which only appear for 4 or more jets, are not included, but 
are expected to contribute a negligible rate (they are fully included in the NJETS code by 
Berends et al [5]). In both initial and final states we only assume as quark types u, d, s 
and c. For the generation of events with heavier quarks (b and t), we suggest using the 2Q 
element of the package. 

As a default, the code generates kinematical configurations defined by cuts applied to 
the following variables: 
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jproc subprocess jproc subprocess jproc subprocess 

1 99 ^ 99 2 gg ^ gg 3 gq -> qg 

4 qg ^ qg 5 gg ^ qq 6 qq ^ qq 

7 qq' — ► gg' 8 gg' — > gg' 9 gg — ► gg 

Table 25: Subprocesses included in the TVjets code. N — 2 additional final-state 
gluons are not explicitly shown here but are included in the code, with TV up to 
6. For example, the subprocess jproc=4, in the case of TV = 4 corresponds to the 
final state qg — > qggg. The details can be found in the subroutine self lav of the 
file Njetslib/Njets.f . For each process, the charge-conjugates ones are always 
understood. 



N jets 


N = 2 


iV = 3 


N = 4 


N = 5 


N = 6 


LHC (nb) 


375(1)-10 3 


24.5(2)-10 3 


4,174(8) 


709(1) 


126.3(4) 


FNAL (nb) 


23,706(60) 


857(5) 


90.9(1) 


8.66(1) 


826(1)-10~ 3 



Table 26: <r(N jets) at the Tevatron and at the LHC, with parameters and cuts 
given in cqs. (3.25)-(3.26). 



The respective threshold values can be provided by the user at run time. Additional cuts 
can be supplied by the user in an appropriate routine. In the code initialization phase, 
the user can select between 2 choices for the parameterization of the factorization and 
renormalization scale Q. A real input parameter (qfac) allows to vary the overall scale 
of Q, Q = qfac x Qq, while the preferred functional form for Qq is selected through an 
integer input parameter (iqopt=0,l). In more detail: 



iqopt 





1 


Ql 


s 


(pt) 



where {pt) is the average pt of the final state jets. 

We cross-checked our calculation against the results for pp — ► iVg (with N < 6) given 
in[52]. Some benchmark results are given in Table 26, obtained for the following set of cuts 
and scale choice: 

Q 2 = (Pt) (3.25) 
p r jct >20GeV, \r}j\ < 2.5, ARjj > 0.7 . (3.26) 



4. Conclusions 

We presented in this paper a new MC tool for the generation of complex, high-multiplicity 
hard final states in hadronic collisions. To the best of our knowledge, a large fraction of 
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the processes we discussed have never been calculated before in the literature to the level 
of jet multiplicities considered here, due to the complexity of the matrix elements involved. 
In addition to the evaluation of the matrix elements, and the possibility of performing 
complete parton-level simulations, the code we developed offers the possibility to carry out 
the shower evolution and hadronization of the partonic final states. In the current version 
we implemented the Les Houches format for the event representation, and developed the 
relative interface with HERWIG. In the future other hard processes (for example including 
emission of real, hard photons) will be added to the list of available reactions and Higgs 
decay to two or four fermions will be included. 

Our code will allow complete and accurate studies of the SM backgrounds to a large 
fraction of the most interesting new physics phenomena accessible at the Tevatron, at the 
LHC, and at future high-energy hadron colliders. 
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A. The contents of the code package 

The code is written in Fortran77, with the part relative to the matrix element evaluation 
available as well in Fortran 90 (see Appendix C). The code package, contained in the com- 
pressed file alpgen . tar . gz, can be obtained from the URL http : //home . cern . ch/mlm/alpgen. 
Unpacking the zipped tarred file alpgen. tar .gz with the command: 

> tar -zxvf alpgen. tar .gz 

will create the following directory structure: 



2Qlib/ 



Njetswork/ herlib/ wcjetwork/ zjetlib/ 



2Qwork/ QQhlib/ 

4Qlib/ QQhwork/ 

4Qwork/ VF90/ 

Njetslib/ alplib/ 



pylib/ wjetlib/ zjetwork/ 

vbjetslib/ wjetwork/ zqqlib/ 

vbjetswork/ wqqlib/ zqqwork/ 
wcjetlib/ wqqwork/ 



The directories labeled 'lib' contain the source codes for the respective processes. 
The user is not supposed to touch them. The code elements which the user will need to 
access and possibly modify to run his own analyses are contained in the directories labeled 
' work ' . More in detail: 
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• The directory alplib/ contains the parts of code which are generic to the evaluation 
of matrix elements using the ALPHA algorithm. The user should treat this directory 
as a black box. When new processes will be added in the future, this part of the code 
should not change. More in detail: 

— alplib/alpgen. f : contains the general structure of the code, preparing the 
input for the matrix element calculation, the bookkeeping of the cross-section 
determination, the event generation, etc. 

— alplib/alpgen. inc: include file, with the necessary common blocks. 

— alplib/ Aint . f , Asu3 . f , Acp . f : the set of programmes necessary for the cal- 
culations of the matrix element, done by the ALPHA algorithm. 

— alplib/alppdf . f : contains a collection of structure function parameterizations; 
some of them require at run time input tables, which are provided as part of 
the package, and stored in the subdirectory alplib/pdf dat/. The command 
file alplib/pdf dat/hvqpdf contains the necessary logical links to all PDF data 
tables. As a default, we already provide a logical link to this file in all /*work 
directories. It is sufficient to issue the pdflnk command within the desired 
working subdirectory to create the necessary logical links, and allow the use of 
all available PDFs. 

— alplib/alputi . f : This program unit contains a histogramming package which 
allows to generate topdrawer [35] files with the required distributions. Examples 
of the use of this package are provided in the default user files *work/*usr . f . 
Users who prefer other histogramming packages, such as HBO OK, do not need 
to link to this file. 

• The directories * lib/ (*=wqq, zqq, wcjet, wjet, zjet, vbjets, 2Q, 4Q, QQh, 
Njets) contain the parts of the code specific to the generation of WQQ+ jets, 
ZQQ+ jets, W+ jets, W + c+ jets, Z+ jets, nW + mZ + IH+ jets, QQ+ jets, 
QQQ'Q + jets QQH+ jets and N jets events. The respective include files with the 
necessary process-dependent common blocks are included in these directories. The 
user should treat these directories as black boxes. 

• The directories *work/ (*=wqq, zqq, wcjet, wjet, zjet, vbjets, 2Q, 4Q, QQh, 
Njets) contain the parts of the code which the user is supposed to interact with, in 
order to implement his own analysis cuts, etc. They contain the files *usr.f, where 
sample analysis routines are provided. These files host the routines in which the user 
can select generation mode, generation parameters (e.g. beam energy, PDF sets, 
heavy quark mass, etc.) as well as generation cuts (minimum pt thresholds, etc.). 
Here the user initialises the histograms, writes the analysis routine, and prints out 
the required program output. This is the only part of the code in which the user 
is supposed to operate, editing the analyses files, and producing and running the 
executable (see next Section). The versions provided as a default contain already 
complete running examples, with the respective command files (input) containing 
sets of default settings (see next Section for more details on running the code). 
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• The directory her lib/ contains the parts of code relevant for the shower evolution 
using HERWIG. In addition to the HERWIG source and include files for version 6.4, 
this directory includes the file atoher.f, which is the interface between the parton- 
level matrix elements and HERWIG, and the file hwuser.f, which includes the main 
driver for the running of HERWIG, and the part of the code where the user can input 
the analyses routines. More in detail: 

— hwuser.f: user initialization of the analysis. Includes the standard HERWIG 
initialization, histogram intialization, analysis routines, etc. The calls to new 
routines hwigup and hwupro create the interface with the generated hard events. 

— atoher.f: this file contains all routines necessary to read in the unweighted 
events produced by the hard matrix element generator. The routine hwigup 
downloads the initialization parameters of the hard process (process type, num- 
ber of partons, beam energy and beam type, etc.), and allows the main HERWIG 
initialization. The routine hwupro is called for each event: it reads the event 
kinematics, flavour and colour information from the file of unweighted events, 
and translates the event data to allow the HERWIG processing of the shower. 
This file should be treated by the user as a black box. 

— herwig64.f: the herwig64 source code. 

— HERWIG64. INC, herwig6400 . inc: HERWIG common blocks 

— pdf dummy. f: dummy PDF routine, required by HERWIG unless the CERN li- 
brary PDF sets are used. As a default, the current version runs with the default 
HERWIG PDF set, regardless of the PDF set which was used to generate the 
hard process. We verified that this does not affect the features of the showered 
final state. Nevertheless we plan in the future to enforce the consistency between 
PDF set used in the hard generation and in the shower evolution. 

• A similar directory and file structure is provided for PYTHIA. 
B. Running the code 

To compile the code for the WQQ process 3 , change directory to: 

> cd wqqwork 

A Makefile is provided for compilation. Issue the comand 

> make wqqgen 

and the executable wqqgen will be prepared. It can be run interactively, inputting from 
the keyboard the run parameters requested by the code, or using the default command file 
input, issuing the command 

3 Analogous procedures allow compilation and run of other processes. 
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> wqqgen < input 



Editing the file input allows to change the initialization defaults (e.g. the number of jets, 
the heavy quark masses, the PDF sets, etc.). Notice that altering some of the inputs in the 
input file may influence the sequence of parameters requested from the code; for example, 
depending on whether light jets are requested or not in a given process, the generation 
cuts for these jets will or will not be required in input. The input file should therefore be 
edited in a consistent manner. 

The first input parameter requested is the running mode imode. The three available 
running modes are discussed in detail in the following subsections. 

The codes vbjets, 2Q and QQh include the option idecay, to allow for the decay of 
the generated on shell vector bosons and top quarks, taking into account spin correlations 
between fermion decay products by means of decay matrix elements. More details on how 
to run the codes are discussed in the dedicated subsection. 

B.l imode=0 

The simplest option is imode=0, where events are generated according to the selected cuts, 
a total cross section is evaluated, and the user can use the routine evtana to analyse the 
event and fill histograms with the desired distributions. To facilitate the job of the user, 
we provide a (redundant) array of kinematical variables relative to the event. The array 
is initialised in the routine usrfll contained in the file wqqlib/wqq. f , and is stored in 
the common block usrevt contained in the include file wqq.inc. Examples of variables 
provided include: 

• pin (4 ,2) : momenta of the incoming partons 

• pout (4, 8) : momenta of the outgoing particles (maximum 8 outgoing partons) 

• pjet(4,8) : momenta of the final-state partons (i.e. quarks and gluons) 

• pt j (8) , etaj (8) : transverse momentum and pseudorapidity of the final-state par- 
tons 

• pbott(4) (pbbar(4)), ptb (ptbb) : momentum and transverse momentum of the 
heavy (anti) quark 

• plep(4) (ptlep) , pnu(4) (ptmiss) : (transverse) momentum of the charged lep- 
ton and neutrino 

• drj j (8,8) : AR separation in rj — <\> space among the final-state partons 

• drbj (8) (drbbj(8)) : AR separation in r] — (j) space between the heavy (anti) quark, 
and the final-state partons 

• etc. . . 

Similar sets of variables are provided for the other processes. As an output the user will 
find the following files: ('file' is the label assigned by the user at run-start time): 
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• file. stat: the header of this file contains information on the run: value of the 
input parameters (EW and strong couplings, beam types and energies, PDF set), 
hard process selected and generation cuts. Furthermore, this file reports the results 
of each individual integration cycle, with total cross sections, as well as individual 
contributions from the allowed subprocesses. It gives the maximum weights of the 
various iterations, and the corresponding unweighting efficiencies, and the value of 
the cross-section accumulated over the various iterations, weighted by the respective 
statistical errors. 

• file . top: includes the topdrawer plots of the distributions, if requested. The default 
normalization of the histograms is in pb/bin. 

• cnfg.dat: file required by ALPHA, generated at run time; it is not needed for the 
analysis of the output, and will be recreated anew any time the code runs, so the 
user should not bother about it, and it can be safely deleted at any time. 

• file.mon: produced/updated after each 100K events; it contains information on the 
status of the run, dumped every 100K events. It is useful to monitor the progress of 
the run. In addition to this monitoring tool, the user can choose to perform other 
tests during the run, in order to save partial information, or monitor the evolution of 
the plotted distributions. Every 100K events the program calls the routine monitor, 
contained in the user file wqqwork/wqqusr . f , where the user can select which oper- 
ations to perform. As a default, the provided routine prints out each 1M events the 
topdrawer file with the distributions being histogrammed. In case of crash, the results 
relative to the statistics accumulated up to that point are therefore retrievable. 

• file . grid*: The phase-space is discretised and parameterised by a multi-dimensional 
grid. During the phase-space integration, a record is kept of the rate accumulated 
within each bin of each integration variable. At the end of an integration cycle ("it- 
eration" ) , the total bin- by-bin rates are used to improve the grid sampling efficiency. 
This is achieved by assigning sampling probabilities proportional to the bin integrals 
(we ensure however that 20% of the sampling is uniformly distributed among all bins, 
to avoid artificial biases introduced by runs with limited statistics). A subsequent 
iteration can then benefit from a better sampling. The state of the grid at the end 
of each iteration is saved in the file file.gridl. Since the first few iterations give 
rise to distributions which are likely to be biased by large statistical fluctuations, we 
separate a phase of grid warm-up from a phase in which events will be generated and 
distributions calculated. The user should then specify in the input file the number 
of warm-up iterations, the number of events to be calculated for each iteration, and 
then the number of events that will be used for the final event generation and for 
the analysis. At the end of the event generation, the grid will be saved to the file 
file.grid2. We keep this separate from file.gridl to allow the user to choose 
whether to start a new generation cycle using the grid status at the end of the pre- 
vious warm-up phase, or at the end of the previous generation phase. These choices 
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are made by the user in the input file, selecting the variable igrid to be (to reset 
the grids and start a new grid optimization), 1 (to start the new run using the grid 
obtained during the previous warm-up phase) or 2 (to start the new run using the 
grid optimised at the end of the latest event generation) . Both grid files are saved in 
f ile .grid*-old (*=1,2) at the beginning of each run, to allow recovery of the grid 
information in case of run crashes or mistakes. 

B.2 imode=l 

Running the code with imode=l offers the same functionality than imode=0, but will in 
addition write the weighted events to a file. To limit the size of the file, only events which 
passed a pre-unweighting are saved. The pre-unweighting is based on a maximum weight 
Wtmp which is equal to 1% of the actual maximum weight at the moment of the generation 
of the event: wt mp = w max /100. An event with weight w passing the pre-unweighting is 
then assigned a weight w' = wt mp if w < wt mp , or w' = w if w > wt mp - The weight w' is 
then saved to a file, together with the random number seed which initiated the generation 
of this event, and with the value of x\ (useful to check the sanity of the file when it will 
be read again for the unweighting) . Some statistical information on the run, including the 
total number of generated events, the integral, and the overall maximum weight, are saved 
as well in a separate file. The file with weighted events is to be used for a later unweighting. 
One can easily verify that the pre-unweighting procedure does not introduce any bias in 
the final unweighting. The random number seed will then be sufficient to regenerate the 
full kinematical, flavour and colour information on the event. The size of each event is 
57 bytes. Make sure you have enough disk space to write out the number of events you 
require. In addition to the files listed above for imode=0, as an output the user will find 
the following files: 

• file. par: includes run parameters (e.g. beam energies and types, generation cuts, 
etc), phase-space grids, cross-section and maximum- weight information; 

• file.wgt: for each event we store the two seeds of the random number generation, 
the event weight, and the value of x\ for the event (as a sanity benchmark after the 
kinematics has been reconstructed from the random seeds) 

As a default, the bookkeeping of the weight distribution is kept in the routine *usr.f, 
and the relative data are printed in the topdrawer file file. top. The study of the weight 
distribution can guide the user to a more efficient choice of maximum weight before starting 
the event unweighting. 

B.3 imode=2 

After a run with imode=l, a run with imode=2 will perform the unweighting of the already 
generated events, and will prepare the input file for the processing of the events with 
HERWIG or PYTHIA. The code reads first the phase-space grids used for the weighted- 
event generation and the maximum weight from file. par. The user has the possibility 
to edit file. par and replace the maximum weight with a different value, if he convinces 
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himself that a more efficient unweighting can be obtained, without biasing the sample, 
by selecting a smaller maximum weight 4 . We are working on techniques to perform this 
optimised unweigthing in an automatic way. The code will then scan the file file.wgt, 
containing the weighted events. A comparison of the event weights with the maximum 
weight is made, and the unweighting is performed. The kinematics of each unweighted 
event is reconstructed from the relative random number seed. The colour flow for the 
event is then calculated, and the full event information is written to a new file. This file 
will be the starting point for the generation of the full shower, to be performed using 
HERWIG or PYTHIA. As an output the user will find the following files: 

• f ile_unw. stat: includes cross sections, max weight, etc; 

• f ile_unw.top: While no new events are generated, the analysis routines used when 
running in imode=l are applied to the unweighted events, and the relative distribu- 
tions are evaluated. In this way the user can compare distributions before and after 
unweighting; 

• file.unw: list of unweighted events, including event kinematics, flavour and colour 
structure, and event weight. 

B.4 Decay of top quarks and vector bosons with spin correlations 

The on-shell top quarks generated by 2Q and QQh undergo a fully exclusive decay in three 
fermions weighted with exact matrix element. When running in imode=0, 1 the information 
on top decay product momenta is stored in the matrix idee (4, 3, 2), available in the 
routines usrfll and usreut. The meaning of the entries is as follows: idee (1 : 4, i , j ) is 
the four momentum (p x ,p y ,p z , E) of the i-th decay product of the j-th particle (jtl and 
jtbl are the labels for top and antitop quark respectively); i=l is the label for the b (b) 
quark; i=2,3 are the labels for the fermion and antifermion coming from the W decay 
respectively. When running in imode=2 the user is required (interactively) to select one 
top decay mode among seven options: 

1 = ev e bb + 2 jets, 

2 = [iv^bb + 2 jets, 

3 = Tv T bb + 2 jets, 

4 = Ivibb + 2 jets, (I = e, /i, r) 

5 = luil'vi/bb, {I = e, /i, r) 

6 = bb + 4 jets, 

7 = fully inclusive. 

4 The user must however avoid using a maximum weight smaller than 1% of the true maximum weight, 
because of the threshold used in the pre-unweighting phase. 



-34- 



Non-zero masses of the VF-decay products are introduced by rescaling their momenta, and 
keeping invariant the W 4-momentum. 

The on-shell W and Z gauge bosons generated in vbjets undergo a fully exclusive 
decay. While the decay products of the W can be changed when running in imode=2 (see 
below), given the universal electroweak couplings of the W to fermions, the Z decay options 
must be specified at the very beginning. When running in imode=0 , 1 the information on 
the W, Z decay product momenta is stored in the matrix idee (4,4, maxpar-2) , available in 
the routines usrf 11 and usreut. The meaning of the entries is as follows: idee ( 1 : 4 , i , j ) 
is the four momentum (p x ,p y ,p z , E) of the i-th decay product of the j-th particle; i=l,2 
are the labels for fermion and antifermion respectively. The ordering in j correspond to 
the ordering in which Z and W's are generated. The flavour of the Z decay products is 
stored in the variable zfl (maxpar-2) according to PDG conventions. For every Z boson 
in the final state the user have to select its decay mode in the input file by entering an 
integer string with the decay modes of the individual Z's according to the following table 

1 = uv, (summed over all flavours), 

2 = / + /~, (summed over all flavours), 

3 = qq, (summed over all flavours), 

4 = bb, 

5 = fully inclusive. 

Concerning the decay modes of the W's, they have to be specified when running in imode=2 
as follows (the code asks for this automatically and interactively): 



1 = 


ev e , 


2 = 




3 = 


TV T , 


4 = 


Wi(l = e,/i,r), 


4 = 


qq, 


5 = 


fully inclusive. 



The same final state options for W bosons are available also for wjets and wqq with 
imode=2. Finite fermionic masses of the decay products are introduced by rescaling the 
momenta and preserving the vector boson 4-momentum. 

B.5 Running HERWIG or PYTHIA on the unweighted events 

A HERWIG executable can be obtained starting from the default driver herlib/hwuser . f . 
To compile and link, issue the command: 

> make hwuser 
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from the herlib directory. The resulting executable, hwuser, should then be run in the di- 
rectory containing the unweighted event file. The default code simply generates the shower 
evolution for the given events, and prints to the screen the particle content of the first 
few events (number of printed events set by the variable maxpr). Analysis code should be 
provided by the user, by filling the standard HERWIG routines hwabeg , hwanal , hwaend. 
A log file f ile-her . log documenting the inputs and outputs of the run is produced. An 
analogous interface for running PYTHIA is available. 

C. Portability and Fortran 90 

The code was tested on several platforms, including Linux based PC's, Digital Alpha 
Unix, HP series 9000/700, Sun work stations and MAC-OSX with a g77 compiler. The 
user should however check that the compilation options provided by default in the first few 
lines of the Makefile files, including the choice of Fortran compiler, are consistent with 
what he has available. We shall be happy to receive comments related to the portability 
of the code, and will update the code to improve its usability. 

Toghether with the old Fortran77 (F77) version of the ALPHA code [36, 27], we provide 
a new Fortran90 (F90) version. The evaluation of the ALPHA matrix elements with the 
F90 version is a factor of five to twenty times faster than the F77 version, depending on the 
selected process (the more complex the process, the better the improvement). When the 
overhead of the rest of the code (phase-space, parton densities, etc) is added in, the overall 
performance of the code improves by a factor of two to five (we stress that only the ALPHA 
part of the code is available in F90; users unfamiliar with F90 should not be discouraged 
from using this version, since this component is a black box, and its use is compatible with 
the F77 part of the code which the user has access to. Furthermore, the F90 executables 
will run using the same input files as the F77 versions of the code, and produce the same 
results, to machine precision). 

To link the F90 version of ALPHA it is sufficient to input the choice of F90 compiler in 
the Makefile, and issue the comand make wqqgen90, which will produce the executable 
wqqgen90. 

For user who do not have access to a F90 compiler, we provide one suitable for running 
on PC's with Linux operating systems. To set it up, proceed as follows: 

• go to the ALPGEN home directory; 

• issue the command make ft90V. This will unpack the file ft90V.tar.gz and in- 
stall the Vast/Verydian F90 compiler into the directory F90V. This software was 
distributed freeware for personal use only by Pacific-Sierra Research. Before use, 
you are therefore supposed to agree with the license term contained in the directory 
F90V/. In the same directory the user can find some documentation on the compiler, 
including the list of supported platforms. 

• Move to the desired work directory (e.g. wqqwork); 

• issue the command make wqqgen90V which will produce the executable wqqgen90V. 
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